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Abstract. We consider the transverse structure and stability properties of relativistic jets formed in the course of 
the collapse of a massive progenitor. Our numerical simulations show the presence of a strong shear in the bulk 
velocity of such jets. This shear can be responsible for a very rapid shear-driven instability that arises for any 
velocity profile. This conclusion has been confirmed both by numerical simulations and theoretical analysis. The 
instability leads to rapid fiuctuations of the main hydrodynamical parameters (density, pressure, Lorentz factor, 
etc.). However, the perturbations of the density are effectively decoupled from those of the pressure because 
the beam of the jet is radiation-dominated. The characteristic growth time of instability is much shorter than 
the life time of the jet and, therefore, may lead to a complete turbulent beam. In the course of the non-linear 
evolution, these fiuctuations may yield to internal shocks which can be randomly distributed in the jet. In the case 
that internal shocks in a ultrarelativistic outfiow are responsible for the observed phenomenology of gamma-ray 
bursts, the proposed instability can well account for the short-term variability of gamma-ray light curves down 
to milliseconds. 
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1. Introduction 

Catastrophic stellar events li ke ma ssive stellar core col- 
lapse (see, e.g., van Putten 2001 , Meszaros 2001 and 



(Paszynski 
al 



1989 



references therein) or merging of a binary neutron star 
1996| , Eichler et al. 
1993 ) have been proposed to 
(GRB) 



Mochkovitch et 
explain the energetics 
Nowadays, an increasing 



of gamma-ray bursts 
amount of observational evidence stresses the association, 
at least, of some GRB events with massive progenitors 
(Bulik, Belczynski 
Hanlon et al 



2000 



Zbijewski [1999|; B loom et al. 
Reeves et al. 



1999 



2002). One of the most 



attractive scenarios is the coUapsar model (Woosley 1993, 
MacFadyen & Woosley 1999) which can produce the re- 
quired energy of a GRB even at cosmological distances. 
This model involves the collapse of the central core of a 
massive, evolved star to a newly- formed black hole (BH). 
The progenitor star can be, for instance, a rotating Wolf- 
Rayet star ( |1993D . Hydrodynamic simulations of collap- 
sars have been performed for a 35Mq main-sequence star 
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whose 14M0 helium core collapses to form 2-AMq black 
hole. Provided that the core has a sufficient amount of an- 
gular momentum, a massive geometrically thick accretion 
disc of several tenths of a solar mass can be formed around 
the BH. The BH accretes matter from the disc at rates of 
the order of solar masses per second. The burst is gener- 
ated by a local energy deposition due to the annihilation of 
ly — D coming from the accretion disk or/and due to the re- 
lease of the BH spin energy by means of magnetic fields. In 
either case, the energy is preferentially released along the 
rotation axis, close to the central engine and gives rise to 
a relativistically expanding bubble of radiation and pairs. 
The duration of the burst is given by the hydrodynamic 
time scale for the core of the star and its Helium enve- 
lope to collapse or by the viscous evolution time for the 
accretion disc, whichever is greater. In order to produce a 
long-duration GRB with a complex pulse structure, accre- 
tion is argued to proceed over a period of time comparable 
to the prompt phase of a GRB. The energy released due to 
accretion is sufficient to drive a collimated, baryon-dilute 
fireball that penetrates the outer layers of the star forming 
a relativistic jet (Aloy et al. 2000). 
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Recent observations indicate that the rapid temporal 
decay of several GRB afterglows is consistent with the evo- 
lution of a highly relativistic jet after it slows down and 
spreads laterally rather than with a spherical blast wave 
(Sari, Piran & Halpern 199!;, Halpern et al. 1999, Kulkarni 



at al. 1999, Rhoads 1999). Therefore, formation of rela- 
tivistic jets of baryon-clean material with bulk Lorentz 
factors ^ lO'^ — lO'^ represents a mayor problem in the col- 
lapsar model of GRB. Detailed simulations of formation 
and propagation of a relativistic jet in coUapsars have been 
performed by Aloy et al. ( pOOO ). Assuming an enhanced 
efficiency of energy deposition in polar regions (with a con- 
stant or varying deposition rate), the authors obtained 
an ultrarelativistic jet along the rotation axis, which is 
highly focused (with the half-opening angle ~ 6 — 10°) 
and is capable of penetrating the star. The simulations 
were performed with the multidimensional relativistic hy- 
drodynamic code GENESIS (Aloy et al. 199E ) using a two- 
dimensional spherical grid. A relativistic jet in this model 
forms within a fraction of a second and exhibits the main 
morphological elements of "standard" jets: a terminal bow 
shock, a narrow cocoon, a contact discontinuity, and a hot 
spot. The maximum Lorentz factor is as large as ^ 25 — 30 
at break-out, and can be even greater after the break-out 
(« 44 at the end of simulations). This Lorentz factor is 
only a bit smaller than the critical value ~ 10^ — 10'^ that 



one requires for the fireball model (Cavallo & Rees 1978 
Piran 119991). 



According to current views, the GRB is made either 
as the jet encounters a sufficient amount of mass in cir- 
cumstellar matter or by internal shocks in the jet (Rees & 
Meszaros |l99||99§ Daigne & Mockkovitch |199S| , |2000D . 
Therefore, the properties of jets from coUapsars are of 
crucial importance for understanding the mechanism of 
GRBs. For example, the non-uniformity across the jet can 
influence the structure of internal shocks and, hence, the 
gamma-ray emission from these shocks. A transverse gra- 
dient can also drastically change the temporal decay rate 
of afterglows, making the decay flatter or steeper depend- 
ing o n the transverse structure (Meszaros, Rees & Wijers 
1998). Gamma-ray light curves also show a great diver- 



sity of time dependences ranging from a smooth rise and 
quasi-exponential decay, through light curves with several 
peaks, to variable light curves with many peaks, and sub- 
structure sometimes down to milliseconds. Hydrodynamic 
instabilities arising in the jet can be responsible for the 
fine time structure of GRBs in a row with a tempo- 
ral variability caused by accretion onto the BH and the 
circumstellar interaction. This complex time dependence 
can provide clues for the understanding of the geometry 
and physics of the emitting regions. Apart from this, in- 
stabilities and their associated fluctuating hydrodynamic 
motions can lead to the generation of magnetic fields 
in a highly conductive e^e^ plasma. In the presence of 
turbulent magnetic fields, the electrons can produce a 



chrotron photons can extend the spectrum even into the 
GeV range (Meszaros, Rees & Papathanasiou 1994). 

In the present paper, we consider the stability proper- 
ties of jets which are formed and propagate in coUapsars. 
Large-scale hydrodynamic instabilities can be the reason 
for the observed morphological complexity of "standard" 
jets, and this has motivated many analytical (see, e.g 
Birkinshaw [1984 ^99l\ ^997[ Har dee & 
Zhao et al. |l992l Hanasz & Sol |1996D 



Norman |1988 , 
and numerical 



(Hardee et al.Jl992| Hardee et al |1998| , Bodo et al. |1998 , 
Micono et al. 2000, Agudo et al. |2001 ) studies of the sta- 
biUty properties of jets. Usually, the jet is considered as 
a beam of gas with one bulk velocity and constant den- 
sity surrounded by a very narrow shear layer separating 
it from the external medium. One possible mechanism 
of destabilizing such jets is often attributed to the well- 
known Kelvin-Helmholtz instability which, in its classical 
formulation, is the instability of a tangential discontinu- 
ity between two flows, generally having different density 
(see, e.g.. Landau & Lifshitz 1978 , Chandrasekhar 1981). 
However, jets in coUapsars show a complex structure with 
the presence of strong transverse shear and substantial 
density stratification (Aloy et al. 2000). Obviously, the sta- 
bility properties of such sheared, stratified jets may well 
be different from those of jets with constant bulk velocity 
and density. 

The outline of this paper is as follows. In Section 2, we 
represent the results of numerical calculations and discuss 
the transverse structure of a jet originating in a coUap- 
sar. In Section 3, we represent the linear stability analy- 
sis of a sheared, stratified jet by making use of a WKB- 
approximation. Finally, in Section 4, we discuss the pos- 
sible role of the instability in the evolution of jets. 

2. Numerical simulations of jets from coUapsars 

The s imulations performed are the same as in Aloy et al. 
( ^000| ). Here we reproduce the technical details for com- 
pleteness. The computational mesh is a two-dimensional 
spherical grid (with coordinates r,9). In the r-direction 
the computational grid has 200 zones spaced logarithmi- 
cally between the inner boundary and the surface of the 
helium star at i? = 2.98 x 10^" cm. Equatorial symme- 
try is assumed and the working angular resolution is 0.5° 
close to the polar region (0° < 9 < 30°) and decreases 
logarithmically between 30° < 6* < 90°. 

A central Schwarzschild BH of mass 3.762Mo pro- 
vides the gravitational field that couples the system. Self- 
gravity of the star is neglected, i.e., we consider only the 
gravitational potential of the BH. Our equation of state 
(EoS) includes the contributions of non-relativistic nucle- 
ons treated as a mixture of Boltzmann gases, radiation, 
and an approximate correction due to e"'"e~-pairs as de- 



synchrotron radiation spectrum (Meszaros & Rees 1993, 
Rees & Meszaros 1994) similar to that observed (Band et 
al. 1993). The inverse Compton scattering of these syn- 



scribed in Witti, Janka & Takahashi (1994). Complete ion- 
ization is assumed, and the effects due to degeneracy are 
neglected. We advect nine non-reacting nuclear species 
which are present in the initial model: G^^, O^^, Ne^°, 
Mg^'^, Si^*, Ni^^, He**, neutrons and protons. 
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In a consistent collapsar model a jet will be launched 
by any process which gives rise to a local deposition 
of energy and/or momentum, as e.g.,z/i/-annihilation, or 
magneto-hydrodynamic processes. We reproduce such a 
process by depositing energy at a prescribed rate homo- 
geneously within a 30° cone around the rotation axis. In 
the radial direction the deposition region extends from 
200 km to 600 km. We have investigated a constant en- 
ergy deposition rate E = 10^° ergs~^, of _E = 10^^ ergs^^, 
and a varying deposition rate with a mean value of 
XQ50g^.gg-i rpj^g constant rates bracket the expected de- 
position rates of collapsar models, while the varying rate 
mimics, e.g., time-dependent mass accretion rates result- 
ing in time-dependent z^i/-annihilation (MacFadyen & 
Woosley |l999| ). 

We have endowed the star with a Gaussian atmo- 
sphere, which at _Ra — 1-8 i? passes over into an external 
uniform medium with a density 10~^gcm~'^ and a pres- 
sure 10~^p{R). The computational domain is extended to 
Rt = 2.5AR with 70 additional zones. At the end of the 
simulation, a relativistic jet has broken out from the sur- 
face of the star and propagates through the added atmo- 
sphere. 

We will focus in the following on the model C50 of 
Aloy et al. (200C). The logarithm of the rest-mass density 
and the Lorentz factor after an evolution time of 5.240 s 
are displayed in Figs. |^ and 0. Note that we have enlarged 
the horizontal scale in Fig. 0in order to improve the read- 
ability of the transverse structure of the jet. 

In Fig. H, we plot the radial dependence of a Lorentz 
factor, r, for a few selected polar angles. This dependence 
shows quite well that the jet is highly coUimated. The 
Lorentz factor reaches maximum values at the axis where 
it can be as large as ~ 20 — 30. However, F drops very 
rapidly with increasing polar angle and, at 9 — 10°, it 
becomes close to 1 in the main fraction of a collapsar vol- 
ume (r > 10® cm). Within the collapsar, on average, the 
Lorentz factor increases with r because of the decreasing 
density. The maximum in F is caused by a very strong rar- 
efaction wave behind the recoUimation shock. Note that 
the Lorentz factor shows a well developed small-scale 
structure in the region with r > lO^cm. The character- 
istic length scale of variations increases from ~ lO^cm at 
r — lO^cm to ~ lO^cm at r = lO^'^cm. These fluctuations 
may indicate the presence of a very rapid hydrodynamic 
instability in the jet. 

Fig. ^ shows the radial dependence of the rest-mass 
density for a few selected polar angles. The density shows 
an even stronger dependence on 9 than does the Lorentz 
factor. For instance, the density contrast between the axis 
(where the density is extremely low, ^ 10~^ — 10~''g/cm^) 
and 9 = 10° is 4-6 orders of magnitude. This is in large 
contrast to the relatively small average change of the den- 
sity with r which does not exceed one order of magnitude 
when r varies from 2 x lO^cm to R. The main reason for 
this change is the fact that the radial grid stretches log- 
arithmically and, therefore, as r grows, there is a loss of 
numerical resolution that may suppress small-scale struc- 



tures. Like the Lorentz factor, the density shows a very 
remarkable small-scale structure approximately with the 
same length scale as F. However, fluctuations of the den- 
sity are larger (up to a factor of 10^ — lO'^). The reason 
for this behavior of the density is that the beam is radi- 
ation dominated according to our EoS (see Sect. ||) be- 
cause the average temperature of the beam is « 5 x 10* K. 
In this radiation-dominated flow, the relative variation of 
the density, resulting from variations of temperature and 
of pressure, can be written as 



Ap 



p pKT p 



(1) 



where T, p and K are the temperature, the pressure and 
the Boltzmann constant, respectively. As the relative vari- 
ations of the pressure and of the temperature are similar, 
it turns out that 



Ap p Ap 
p pKT p 



(2) 



Therefore, p/{pKT) >> 1, we conclude that the rela- 
tive variations of the density must be much larger than 
the relative variations of pressure or of temperature in a 
radiation-dominated beam. 

In Fig. we show transverse dependences of the radial 
velocity, sound speed, density and Lorentz factor at differ- 
ent distances. All quantities are coUimated in a cone with 
a half-opening angle of about 6-10° with a better coUima- 
tion for the Lorentz factor and a worse one for the density. 
Only close to the black hole (r ~ 10* cm) is the colli- 
mation of the represented quantities (except the Lorentz 
factor) worse, ~ 20°. Note also that all quantities except 
the Lorentz factor exhibit almost a self-similar behavior 
at r > 10*cm. The velocity profiles show the presence of a 
strong shear which is crucial for the dynamics of jets and 
their stability properties. The sound speed is close to the 
maximum possible value, c/VS, in a significant fraction 
of the jet volume, i.e., the jet plasma is relativistic in the 
sense of both the bulk velocity and the thermal content. 
The mean temperature is ~ 5 x 10*K, and the pressure is 
dominated by radiation. That is why the sound speed is 
relativistic despite the temperature being non-relativistic 
for particles. Note the presence of a substantial transverse 
gradient of the sound speed that may play an important 
role in the formation of internal shocks. Jets in coUapsars 
exhibit a highly inhomogeneous density in contrast to the 
standard jet models. The density becomes extremely low 
at the jet axis and it grows almost exponentially with the 
angle reaching larger values at the jet surface. This behav- 
ior is caused by a velocity profile decreasing with 9 and, 
likely, is typical for all sheared jets. 

3. The dispersion equation for perturbations 

We model the jet from coUapsars by an infinitely long 
plasma cone with a half-opening angle ^o- From simula- 
tions, the typical value of 9o is somewhat around 6-10°. 
Plasma inside the jet moves with a velocity V = V{r, 9)er 
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with respect to the ambient medium; r, ^, tp are the spher- 
ical coordinates with e^, e^, being the corresponding 
unit vectors. As seen from the behavior of the velocity 
profile (Fig. ||), we can distinguish two different regions 
within the jet: the central core where the flow is ultrarel- 
ativistic, and the surrounding transition layer where the 
velocity decreases from relativistic to non-relativistic val- 
ues. In the core region, the Lorentz factor is large, > 3, and 
may even reach values as high as 25 — 30 (see Fig. ^. 
It is convenient to distinguish between these two regions 
in our analysis because of the difficulty of describing both 
relativistic and non-relativistic flows within the framework 
of the same analytical formalism. Basically, the thickness 
of the transition layer is smaller than the radius of the 
core. Since our main interest is the processes in the core 
region we will model the effect of the transition layer in 
terms of the boundary conditions. Note that outside this 
transition layer there exists an extended region where the 
flow is non-relativistic but we neglect the influence of this 
region because of its low energy. 

In relativistic hydrodynamics the continuity, momen- 



tum, and energy equations read (see, e.g., Weinberg 1972) 



dp 
dt 



'di 



V • (Fv) 



'di 



+ (v • V)v 



V dp 



0, (3) 



(4) 



-(pn-^)+vV(pn-^) = 0, 



(5) 



where F = (1 — w^/c^)"^/^ is the Lorentz factor, p and 
p are the gas pressure and density, respectively; n is the 
number density in the fluid's rest frame; 7 is the adiabatic 
index. 

Our analysis of stability is based on the linearized set 
of equations Small perturbations will be marked 

by the index 1; for unperturbed quantities subscripts will 
be omitted with the exception of vector components. We 
consider the instability which arises on a time scale much 
shorter than the characteristic evolution time scale of the 
jet. Therefore, we can adopt a quasi-steady approximation 
neglecting the time dependence of unperturbed quantities. 
In this approximation, the linearized continuity, momen- 
tum and energy equations are 

F (pi -f vi • Vp + V • Vpi) + [V • (Fvi)-f 



Vlr + V 



dvir 

dr 



d 



(6) 



F>* [vi (vi • V)V + (V • V)vi] + Vpi + -^px 
dp ( pi 2T^Vp, 

o~ Pi + ^ "I 9 ""ir 

p* dr \ & & 



Pi 



VVpi + ^ 



dp 
dr 



V • (Fvi) + 



Vlr 



V 



dvir 

dr 



-Vlr 



^pV d 
dr 



p^^lir'yn (8) 



where p* = p + p/c^. In equations (||)-(^) we took into 
account that Fi = F^(V-vi)/c2. Since a quasi-steady state 
approximation is adopted, we have for the unperturbed 
state F2p*(V • V)V = — Vp, and hence only the radial 
component of Vp is non- vanishing. 

In a conical jet with a small opening angle, the charac- 
teristic radial length scale of unperturbed quantities is of 
the order of the spherical radius and is much greater than 
the transverse length scale. Therefore, we can treat the 
radial and ^-dependences of the perturbations separately. 
If we consider perturbations with the radial wave vector, 
fc, satisfying the condition 



kr > 1, 



(9) 



then the terms on the right hand sides of equations (gl)-(§l), 
which contain only the radial derivative of unperturbed 
quantities, are smaller by a factor kr than the correspond- 
ing terms on the left hand sides containing the derivatives 
of perturbations. Therefore, all terms on the r.h.s. of equa- 
tions (H)-(H) can be neglected if inequality (||) is fulfilled. 
The assumption (^) corresponds to the so-called local ap- 
proximation in a stability analysis. Under this assumption, 
we have from (^)-(^ 



r(pi 



Vl 



Vp + V • Vpi) 

Vlr 



V 



P* [v(Fvi; 

dvir ' 



dr 



V 



F>, [vi + (vi • V)V + (V • V)vi] + Vpi + —pi 



(10) 



(11) 



Pi + VVpi 



TP 
F 



V(Fvi) 



Vlr + V 



dvir 



= 0. (12) 



Since the unperturbed quantities depend neither on t nor 
on if , the dependence of all perturbations on i, r and 
if can be taken in the form exp(iix't — ikr — irmp) if we 
adopt condition (^. In the present paper, we consider ax- 
isymmetric perturbations with m = 0. Substituting this 
dependence into equation (|l^), we can express the pertur- 
bation of pressure in terms of vi , 



Pi 



i"fp 



V • (Fvi) 



-Vvir 



(13) 



where a — uj — kV. The quantity kV is an advective fre- 
quency and has a kinematic origin. The perturbations of 
the velocity can be expressed in terms of pi from the mo- 
mentum equation (|l^) which under our assumptions can 
be rewritten as 



(7) p,F (iCTVi+e^ — — 



(14) 



Aloy et al.: Stability analysis of jets from coUapsars 



5 



Substituting expressions for vir and vw into equation (|13|) 
and taking into account that for a small opening angle 
sinO w 9, we obtain the equation containing pi alone 

2r^V' ( , u:^ 



P_ 

P* 



c? V 



Pi=0- 



(15) 



Here Cg — \/jp/p* is the sound speed, x — rO, and 
the prime denotes a transverse derivative; for example, 
p'l = dpi/dx where dx — rdO . Equation ( p^ ) represents 
the behavior of small perturbations for any velocity and 
density profiles within a conical jet with a small opening 
angle. In our approach, the radial dependence of unper- 
turbed quantities enters in equation ( |l5|) parametrically. 

For our purposes, it will be convenient to use another 
form of the equation (|l5|). Making the substitution 



Pi = 



exp 



k — ] dx 



equation ( [LSD can be transformed into 
where 



(16) 



(17) 



c \ c 



1 

p'V 

p* 



V 

c 

{xvy 



kc/ aV" , and Q is given by 
2 



c 



{xpj 
2xp^, 



(18) 



Since in the present paper we consider a stratified sheared 
jet the expression ( p7[ ) for d iffers from the correspond- 
ing equation derived by Urpin ( ^002 ) only due to the pres- 
ence of terms containing derivatives of p. 

In the core region where the temperature is high and 
the pressure is determined by radiation, the sound speed, 
Cs, is of the order of and, hence, the acoustic fre- 

quency is comparable to that of the light, ~ kc. Since the 
turnover time scale, characterizing the rate of hy- 

drodynamical processes associated with shear is also very 
short, we can expect that typical frequencies in our model 
are sufficiently large. Therefore, we can try to obtain the 
solution of equation ( p^ ) in a high frequency limit when 
^ < 1, or 



a |> fcc/r^ 



(20) 



We will see that this inequality is well fulfilled in our jet 
model. Restricting ourselves to the lowest terms in a small 
parameter ^, we obtain from the expression (|l8|) 



VT2 

p* 



p'V {xvy 



(V) 

2xp^, 



where rj — \ — c^y^/c^. Since in the ultrarelativistic limit 



we have 



c /3 and V the factor 77 varies in a 



narrow range from ^ 1 to 2/3. 

The behavior of q is determined by the 0-dependence 
of the jet velocity and density. We consider the case closely 
related to the calculations represented in Section 2. These 
simulations indicate that the density profile within the jet 
core can be fitted with a sufficient accuracy by a simple 
exponential dependence, 



p{r) w poe" 



(22) 



where po is the density at the jet axis, and Lp is the char- 
acteristic length scale of stratification. Since the density 
contrast between the surface and axis is of the order of 
3-5 orders of magnitude depending on the distance from 
the center of a coUapsar, we can estimate Lp ~ O.Ixq. 

The velocity profile can be modeled by the parabolic 
dependence. 



V{r)^Vo{l~x^/Ll), 



(23) 



where Vb is the velocity at the jet axis, and L„ is the length 
scale of the velocity distribution, L^, ~ a^o- 

Since is large in the jet core, the contribution of the 
last four terms on the r.h.s. of equation ( pT| ) which are pro- 
portional to F^ or F" is small compared to the second term 
(cx F*) if we assume that F^ 3> max(Li,/ip, L^j ^ xLp^ 
in the main fraction of the core region where instability 
arises. After these simplifications, we obtain the following 
expression for q^ in the jet core 



r]a 



(24) 



Consider the solution of equation ( |17D by making use 
of a WKB-approximation which is well justified in those 
regions of the jet where the condition \qx\ > 1 is fulfilled. 
With the accuracy in linear terms in a small parameter 
|ga;|~^, a WKB-solution of equation ( p7| ) can be written 
as 

1 



i J qdx' 



(25) 



(Landau & Lifshitz 1981) where Ci and C2 are constant 



that have to be chosen in such a way as to satisfy the 
boundary conditions. In some cases, the condition \qx\ ^ 
1 holds everywhere within the jet, and the expression (25) 
applies at xq > a; > 0. However, since the dependence 
of g on X is rather complex, it is generally also possible 
that the condition of applicability can break at some point 
X — X:t < xq where q ~ 0. By analogy with quantum 
mechanics this point can be called the turning point. We 
address this case as the most interesting one from the 
point of view of applications to jets from coUapsars. Note 
that we have at a; = a;* 



r] \ c y 



fV- 



Since generally lu is complex, 
(jj = fl — i/r. 



(26) 



(27) 
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where is the frequency of an unstable mode and r is its 
growth time, the condition (p6|) can be satisfied only if 



l/r> \n-kv\ 



(28) 



at the turning point. We are interested in the most pow- 
erful instability (if it exists) with the growth time shorter 
than For such instability the condition ( p8| ) should 

be fulfilled everywhere within the jet. If equation 
holds then is approximately real and < at 
a; > and > at xo > x > x*. 

It is convenient to represent a WKB-solution 
both sides of the turning point in the form 

1 



rx 



qdx' 



qdx' 



(29) 



and Cj'^'''' are con- 



(Landau & Lifshitz |l98l|) where 
stant; indexes r and I refer to the right (x > a;*) and left 
(a; < a;*) from x = x*, respectively. Note that generally 
&;^^C^^ andC^'^Vcf . 

To obtain the eigenvalues and eigenfunctions of equa- 
tion (|l^) one needs the boundary conditions for pertur- 
bations. One of the boundary conditions is obvious: since 
Pi is a physical quantity it should be finite everywhere in- 
cluding the jet axis and, hence, / should vanish at a; = 0. 
If the turning point is not very close to 0, then usually it 
is sufficient to choose an exponentially decreasing solution 
at X < x*. It is known (see, e.g.. Landau & Lifshitz 1981) 
that the solution ( |29| ) at x > matches an exponentially 
decreasing solution beyond the turning point (at x < x*) 
if Cf ^ = (C/2) exp(i7r/4) and C^'^ = (C/2) exp(-i7r/4) 
where C is constant. Hence, the solution at x > x, satis- 
fying the true boundary condition at x = has the form 

/(x) = cos ( r q{x')dx' + ^] . (30) 



The formulation of a true boundary condition at the 
outer boundary can face some problems but, fortunately, 
the properties of modes with a relatively small length scale 
in the 0-direction are not sensitive to the particular choice 
of boundary conditions. We consider the case when the 
density in a surrounding medium is much larger than that 
in the jet. In this case (see, e.g., Glatzel 1988, Wu & Wang 



1991), a much denser surrounding medium plays the role 



of a wall and the perturbations of the displacement in the 
^-direction vanish at the outer boundary, x = xq. As it 
follows from equation (p^, vig vanishes if dpi/dx = at 
X = xq. Using the expression (p^, we can rewrite this 
condition as 



/' + 



r' / c ^ 



f = o. 



(31) 



2/3* 2x0 

Since we consider the modes satisfying condition (|2^) and, 
in a WKB-approximation, the wavelength of perturba- 
tions is assumed to be much shorter than the characteristic 
length scales of any unperturbed quantity, the boundary 
condition at x = xq can approximately be represented as 



f'ixo) « 0. 



(32) 



With this boundary conditions, the dispersion equation 
reads 



q{x)dx = 7ra, 



(33) 



where a — n ^- 3/4, n is integer. Substituting the expres- 
sion (|2j) for q in (|3^), we have 



2:0 /^„2 y2p2 \ 



dx = 



(34) 



We can estimate the characteristic roots of equation 
( ^4[ ) using the mean value theorem. We have 



an 



(35) 



where Ax = xq — x, and the subscript a means the value 
of a quantity at x = Xq, Xq is a mean point within the 
range from x» to xq. Since under the condition ( p8| ) we 
have (7^ ~ — l/r^, then the growth rate is given by 



1 l/'^r^ / a 



2 _ (Ax_)Hf 



where 



(36) 



(37) 



If n < — 3/4, then is positive and perturbations with 
such n are unstable. On the contrary, perturbations with 
very large n satisfying the condition n > A^ — 3/4 are 
stable. Note that A^ '-^ . In the case n N, we have 



(38) 



At Cs ~ c/ \/3, we obtain a very simple estimate of the 
growth time 



(39) 



Jets with a relativistic sound speed turn out to be unstable 
for any velocity profile. The growth time is typically very 
short and depends on the shear being shorter for larger 
\V'\. Also, the instability arises faster for jets with a very 
high Lorentz factor. Note, however, that the origin of the 
factor Ta in equation ( ^ ) is purely kinematic. The growth 
time of the instability in the comoving frame is of the 
order of \V'\, but since we consider the instability in the 
rest frame, the Lorentz transformation of the "comoving" 
growth time yields the result (^) . 

The characteristic growth length, L = Vt, can be es- 
timated, for Cs ^ as 



TalK'l r,' 



(40) 



and is much shorter than the radius of a coUapsar. 
Therefore, the instability can manifest itself even at the 
very early stage of the evolution of the jet. 
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4. Discussion 

We have treated numerically and analytically the instabil- 
ity that can arise in jets from coUapsars. The instability is 
caused by a combined action of shear, which is unavoidable 
in such jets, and an extremely high compressibility asso- 
ciated with a relativistic sound speed. It turns out that 
jets from coUapsars are more unstable than, for example, 
standard extragalactic jets because of the relativistic com- 
pressibility. The fact of instability itself docs not depend 
on a particular shape of the velocity profile: the instability 
can arise for any dependence V{9). However, the growth 
time is shorter for fiows with a stronger shear. Note that 
only non-homogeneous perturbations in the radial direc- 
tion {k 0) can be unstable. 



Although it is commonly believed that three di- 
mensional studies of the stability of relativistic jets 
will include additional instable modes (in many 
cases even more unstable that the axisym metric 
ones) 



199S 



the recent v^rork of Hardee &: Rosen 
has pointed the fact that shear leads both to an en- 
hancement of the axisymmetric modes and a sup- 
pression of the asymmetric modes. Hence, a fully 
three dimensional study would not yield to radi- 
cally different results from the ones that we obtain 
here. 



manifest itself as that where the condition tq ^ r/c is 
fulfilled. Substituting tq, we obtain that 



(5-10)-^^lC 



(41) 



at this radius. At small r {^-^ 8 x 10® cm well below jet 
break-out), the jet is less coUim ated (^o ~ 10 — 12°) than 
on average (see Aloy et al. 2000| ). Therefore, the instability 
should manifest already in regions where F ^ 1 — 2. This 
conclusion is in qualitative agreement with what is shown 
in Fig. 1^ where appreciable fluctuations appear in regions 
with a Lorentz factor between 1 and 2. 

Note that the radial wave vector of unstable perturba- 
tions, k, should satisfy rather restrictive conditions. First, 
k has to be sufficiently large for the applicability of the lo- 
cal approximation (see equation (|^)). On the other hand, 
condition ( p8| ) which is necessary for the existence of the 
turning point, implies that 



V 



xo 



(42) 



If the inequality (^2|) holds then the condition ( ]20| ) used 
in our calculations is certainly fulfilled and the parameter 
^ is small. Combining equations (||) and (^2|) and taking 
into account that xq = tOq, we obtain that 



r, 



- > fcr > 1 



(43) 



In the main fraction of the jet volume, the instability 
grows very rapidly. If we estimate the average Lorentz fac- 
tor as ~ 10 and the half-opening angle as ~ 5° then the 
growth time ( p9| ) is of the order O.Olr/c. This is much 
shorter than the life time of the jet even at r ~ i?. 
Therefore, we can expect that shortly after jet formation 
the instability will generate well developed turbulent mo- 
tions with substantial fluctuations of the Lorentz factor, 
density, pressure, etc. This conclusion is in good agree- 
ment with the results of numerical simulations (see Figs. ^ 
and^). During the jet's propagation, fluctuations can be- 
come noticeable at a relatively early evolutionary stage 
and, hence, at a small distance from the formation region 
because the growth time is sufficiently short even in the 
inner region of the coUapsar. Of course, the initial ampli- 
tude of fluctuations in the jet is quite uncertain. In our 
case, the initial fluctuations are triggered by numerical 
reasons, but it is quite likely that they mimic the irregu- 
larities of the process of accretion onto the central black 
hole. Usually the time tq required for an instability to am- 
plify the amplitude of perturbations to a noticeable value 
is longer than the growth time, r, by some factor which is 
typically ^ (5 — 10), so ro ~ (5 — 10)t. After this time, the 
fluctuations have grown by a factor ^ 10^ — 10^ compared 
to their initial value which probably is sufficient to become 
significant. In the region where tq is shorter than the prop- 
agation time scale, which can be estimated as ^ r/c, fluc- 
tuations reach noticeable values. On the contrary, fluctu- 
ations seem to be insigniflcant in regions where tq > r/c. 
We can define the radius where the instability starts to 



for unstable perturbations. Estimating F '--^ 10 and 
5 — 6° in the main fraction of the jet volume, wc have 
100 '3> kr ^ 1. Numerical simulations (see Fig. ||) indicate 
that the characteristic length scale of fluctuations depends 
on r and varies within the range from O.lr to 0.5r in good 
agreement with the prediction of the theory. Note that 
fluctuations with a shorter length scale cannot be resolved 
with the computational grid used in the simulations. 

All this allows one to speculate that the calculated fluc- 
tuations of parameters within the jet are physical (i.e., 
not simply numerical artifacts) and reflect the presence 
of a very strong shear-driven instability. We can expect 
that inhomogeneities caused by this instability will pro- 
duce shocks in the course of their non-linear evolution 
when faster fluctuations try to overtake slower ones or 
when fluctuations moving in the positive and negative 
radial directions collide. If this is the case then internal 
shocks might be more or less randomly distributed and 
oriented within the jet forming filamentary structures. It 
is often supposed that shocks in a ultrarelativistic wind or 
jet are responsible for GRBs themselves whereas the im- 
pact against the ambient matter of this wind produces an 
external shock which likely produces the observed after- 
glows (Rees & Meszaros 1992, 1994). Shocks can convert 



a portion of kinetic energy into a non-thermal gamma/X- 
ray transient emission which is usually ascribed to par- 
ticle acceleration by shocks. Typically, the efficiency of 
this conversion is not high, ~ 1 — 2 %, but it can be 
much greater (~ 20 — 40 %) for the interaction of fluc- 
tuations with very different Lorentz factors (Kobayashi, 
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Piran & Sari 1997, Kobayashi & Sari 2001). Since in our 
model the calculated fluctuations move with substantially 
different Lorentz factors we can expect a highly efficient 
transformation of their kinetic energy into radiation. The 
proposed instability can also accounts for the rapid vari- 
ability of the gamma-ray light curves, which lasting from 
tens to hundreds of seconds, exhibit variability sometimes 



down to milliseconds (Fishman & Meegan 1995). Likely, 
the most rapid temporal variability associated with the 
shear-driven instability has a time scale of the order of 
the growth time (p9[). At the surface of a coUapsar, for 
example, this time scale is as short as ~ 10~^s. Since the 
jet is highly inhomogeneous and the Lorentz factor varies 
strongly during the jet's propagation, a slower variability 
could also be represented in the gamma-ray light curves. 

Another remarkable inference from the considered 
model is that the turbulent motions caused by the in- 
stability may also be important for the electron-proton 
energy exchange and, particularly, for the generation of 
the magnetic field in jets from collapsars. 
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Fig. 1. Snapshot of the logarithm of the rest-mass density after an evolution time of 5.24 s. The color scale indicates 
the value of the logarithm of the rest-mass density (in CGS units) . The solid arc marks the limit of the exponentially 
decaying atmosphere while the dashed one marks the surface of the collapsing star. The numbers around the box 
represent distances in units of 10^° cm. 
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Fig. 2. Snapshot of the Lorentz factor after an evolution time of 5.24 s. The color scale indicates the value of the 
Lorentz factor at every point of our computational domain. The solid arc marks the limit of the exponentially decaying 
atmosphere while the dashed one marks the surface of the collapsing star. The number around the box represent 
distances in units of 10"'^'' cm. Note that the horizontal scale has been enlarged by a factor of 3 in order to display more 
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Fig. 3. Lorentz factor as a function of the radius for three different angles {9 « 0°, 5°, 10°; see legends in the figure) 
after an evolution time of 5.24 s. 
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Fig. 4. Rest-mass density as a funetion of the radius for three different angles {9 » 0°,5°,10°; see legends in the 
figure) after an evolution time of 5.24 s. 
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Fig. 5. Plots of four physical quantities as a function of the polar angle, 9, at four different distances (10^, 10^, 3 x 10^ 
and 10^° cm; see legends in the figure). The upper left panel corresponds to the radial velocity (in units of c); the lower 
left panel to the rest mass density and, the upper and lower right panels correspond to the sound speed (in units of 
c) and to the Lorentz factor, respectively. The evolution time is 5.24 s. 



